Subtlety of Determining the Critical Exponent u of the Spin-1/2 Heisenberg Model 
with a Spatially Staggered Anisotropy on the Honeycomb Lattice 
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Puzzled by the indication of a new critical theory for the spin-1/2 Heisenberg model with a 
spatially staggered anisotropy on the square lattice as suggested in [3| , we study a similar anisotropic 
spin-1/2 Heisenberg model on the honeycomb lattice. The critical point where the phase transition 
occurs due to the dimerization as well as the critical exponent v are analyzed in great detail. 
Remarkly, using most of the available data points in conjunction with the expected finite-size scaling 
ansatz with a sub-leading correction indeed leads to a consistent v — 0.691(2) with that calculated 
in [l[. However by using the data with large number of spins A'^, we obtain f = 0.707(6) which 
agrees with the most accurate Monte Carlo 0(3) value v = 0.7112(5) as well. 

PACS numbers: 



I. INTRODUCTION 



The discovery of high temperature (high Tc) super- 
conductivity in the cuprate materials has triggered vig- 
orous research investigation on spin-1/2 Hcisenberg-type 
models, which are argued and believed to be the correct 
models for the undoped precursors of high Tc cuprates 
(undopcd antifcrromagncts). In addition to analytic re- 
sults, highly accurate first principles Monte Carlo studies 
of Hcisenberg-type models are available due to the fact 
that these models on geometrically non-frustrated lat- 
tices do not suffer from a sign problem, which in turn 
allows one to design efficient cluster algorithms to sim- 
ulate these systems. Indeed, thanks to the advance in 
numerical algorithms as well as the increasing power 
of computing resources, the undoped antiferromagnets 
are among the quantitatively best-understood condensed 
matter systems. For example, the low-energy parame- 
ters of the spin-1/2 Heisenberg model on the square lat- 
tice, which are obtained from the combination of Monte 
Carlo calculation and the corresponding low-energy ef- 
fective field theory, are in quantitative agreement with 
the experimental results [1]. 

Anisotropic Heisenberg models have been studied in- 
tensely during the last twenty years. Besides their phe- 
nomcnological importance, they are of great interest from 
a theoretical perspective as well [1, 0, [E H ■ ^or ex- 
ample, due to the newly discovered pinning effects of the 
electronic liquid crystal in the underdoped cuprate super- 
conductor YBa2Cu3 06.45 H, 0, the Heisenberg model 
with spatially anisotropic cou pling s Ji and J2 has at- 
tracted theoretical interest [lO|, Further, numerical 
evidence indicates that the anisotropic Heisenberg model 
with staggered arrangement of the antiferromagnetic cou- 
plings may belong to a new universality class, in contra- 
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diction to the theoretical 0(3) universality predictions 
In particular, the numerical Monte Carlo studies 
of other anisotropic spin-1/2 Heisenberg models on the 
square lattice, for example, the ladder and the plaquctte 
models, seem to be described well by the 3-dimensional 
classical Heisenberg universality class [13, even at 
small volumes. These studies indicate that the staggered 
arrangement of the anisotropy in the spin-1/2 Heisenberg 
model might be responsible for the unexpected results 
found in [l|. In order to clarify this issue further, based 
on the expectation that the honeycomb lattice is one of 
the candidates to investigate whether such a new univer- 
sality class does exist, we have simulated the spin-1/2 
Heisenberg model with spatially staggered anisotropy on 
the honeycomb lattice. 

The motivation of our study is to investigate^ whether 
the unconventional critical behavior found in |lj] can be 
observed again by considering the Heisenberg model with 
a similar anisotropic pattern on the honeycomb lattice. 
In particular, one can expect such kind of investigation 
might shed some light on understanding the discrepancy 
between v = 0.689(5) determined in [l| and the most ac- 
curate Monte Carlo value i' = 0.7112(5) for the 0(3) uni- 
versality class [l3|. To achieve this goal, the critical be- 
havior of the spin-1/2 Heisenberg model with a spatially 
staggered anisotropy has been investigated in great detail 
in this study. In particular, the transition point driven 
by the anisotropy as well as the critical exponent ly are 
determined with high precision by fitting the numerical 
data points to their predicted critical behavior near the 
transition. We focus on the critical exponent 1^ so that we 
can perform a more thorough analysis and obtain data 
with as large number of spins as possible. Although we 
find a consistent u = 0.691(2) with that calculated in 
[l[ by using most of the available data points as well as 
taking into account the correction term in the finite-size 
scaling ansatz, we obtain v = 0.707(6) which agrees with 
the known 0(3) Monte Carlo value as well when only the 
data points with a large number of spins are used in the 
fits. 

This paper is organized as follows. In section [Hi the 



FIG. 1; The anisotropic spin-1/2 Heisenberg model investi- 
gated in this study. 



anisotropic Heisenberg model and the relevant observ- 
ables studied in this work are briefly described. Section 
mil contains our numerical results. In particular, the cor- 
responding critical point as well as the critical exponent 
V are determined by fitting the numerical data to their 
predicted critical behavior near the transition. Finally, 
we conclude our study in section HVl 



II. MICROSCOPIC MODEL AND 
CORRESPONDING OBSERVABLES 

In this section we introduce the Hamiltonian of the 
microscopic Heisenberg model as well as some relevant 
observables. The Heisenberg model considered in this 
study is defined by the Hamilton operator 



H — J Sx ■ Sy + J' Sx' ■ Syl , 



(1) 



{xy) 



{x'y') 



where J' and J are antiferromagnetic exchange couplings 
connecting nearest neighbor spins {xy) and {x'y'), re- 
spectively. Figure 1 illustrates the Heisenberg model de- 
scribed by eq. ([1]) on the honeycomb lattice with periodic 
spatial boundary conditions implemented in our simula- 
tions. The dashed rectangle in figure 1, which contains 4 
spins, is the elementary cell for building a periodic hon- 
eycomb lattice covering a rectangular area. For instance, 
the honeycomb lattice shown in figure 1 contains 3x3 
elementary cells. The honeycomb lattice is not a Bra- 
vais lattice. Instead it consists of two triangular Bravais 
sub-lattices A and B (depicted by solid and open circles 
in figure 1). As a consequence, the momentum space of 
the honeycomb lattice is a doubly-covered Brillouin zone 
of the two triangular sub-lattices. To study the criti- 
cal behavior of this anisotropic Heisenberg model near 
the transition driven by the dimerization, in particular 
to determine the critical exponent v, several relevant ob- 
servables are measured in our simulations. A physical 
quantity of central interest is the staggered susceptibility 
(corresponding to the third component of the staggered 




FIG. 2: The space-time volume dependence of the staggered 
susceptibilities Xs of the Heisenberg model in this work for 
various values of anisotropy J' / J. Here stands for the 
number of spins. The lines are added to guide the eye. 



magnetization M^) which is given by 



L1L2 Jq 

= tV dt ^TT[MU0)M!{t)cxp{-(3H)].{2) 
Jq ^ 

Here (i is the inverse temperature, L\ and L2 are the 
spatial box sizes in the 1- and 2-direction, respectively, 
and Z = Tr exp(— /3i/ ) is the partition function. The 
staggered magnetization order parameter Ms is defined 
as Ms = Y.x{-^Y^x- Here (-1)^ = 1 on the A sublat- 
tice and (— 1)"^ = — 1 on the B sublattice, respectively. 
Other relevant quantities calculated from our simulations 
are the spin stiffnesses in the 1- and 2-directions 



(3) 



here i G {1,2} refers to the spatial directions and is 
the winding number squared in the i-direction. Finally, 
the second-order Binder cumulant for the staggered mag- 
netization, which is defined as 



O2 



(4) 



are measured as well. By carefully investigating the spa- 
tial volume (and possibly the space-time volume) and the 
J'/ J dependence of these observables, one can determine 
the transition point as well as the critical exponent v with 
high precision. 



III. DETERMINATION OF THE CRITICAL 
POINT AND THE CRITICAL EXPONENT v 

The simulations done in this work were performed us- 
ing the loop algorithm in ALPS library 17[. In order 
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FIG. 3: The J'/J dependence of psiL (a) and Q2 (6). The hnes are added to guide the eye. 



to estimate where in the parameter space J'/J the tran- 
sition occurs, we investigate the space-time volume de- 
pendence of the staggered susceptibility Xs- To be more 
precise, while Xs grows with increasing space-time vol- 
ume in the broken phase, in the symmetric phase, Xs 
should saturate for large space-time volumes. We fur- 
ther scale P with the square root of the number of spins 
N in order to use this criterion of space-time volume de- 
pendence of Xs- Using this criterion, we obtained figure 
[21 The trend of the curves for different ratios J'/J in 
figured] indicates that the transition takes place between 
J'/J = 1.73 and J'/J = 1.74. Comparing this observa- 
tion with the related critical point calculated in [l| for 
the square lattice case, one sees that the honeycomb lat- 
tice has a smaller value of {J'/J)c. This is reasonable 
since the antiferromagnet on the honeycomb lattice has 
a smaller value of the staggered magnetization density 
Aks per spin than that on the square lattice 0, [l^. As 
a result, it is easier to destroy the antiferromagnism on 
the honeycomb lattice than on the square lattice by in- 
creasing the anisotropy J'/J. Therefore one expects a 
smaller critical value of {J'/J)c on the honeycomb lat- 
tice. Further, since phase transition only occurs for a 
system with infinite degrees of freedom, to investigate 
the critical behavior of a system on a finite volume, it is 
useful to employ the idea of finite-size scaling (FSS) hy- 
pothesis for certain observables. Pioneering by Fisher in 
the seventies and generalizing later by Fisher, Barber and 
others including the derivation of FSS hypothesis from 
the framework of renormalization group 0, HO, HH , 
the FSS hypothesis has become the essential method to 
study the thermodynamic property of a finite system near 
the phase transition. The characteristic length for a fi- 
nite system is the correlation length ^. Consequently it 
is natural to consider its thermodynamics property near 
the transition as a function of the variable i/^, where 
L is the physical size of the system. The FSS hypothe- 
sis states that near the transition, any relevant property 
C'i(i) of a finite system with critical exponent k can be 
described by OL(t) = L'^^^ go{L/ S,)^ where t is given by 



t ~ U~jc)/jc with j ~ J'/J, V is the corresponding criti- 
cal exponent for ^, and go is a smooth function. Further, 
since ^ diverges as ^ \i ~ jc\~'^ near jc, one arrives at 
C'L(t) = L'^^^goiLt") = L'^/^g'^itL^^"), where g'^, is an- 
other smooth function. FSS hypothesis can be rigorously 
derived by considering the scaling of the relevant direc- 
tions in the renormalization group fiows. They can be 
understood intuitively as well by assuming the fiuctua- 
tions of any relevant variable should be invariant at all 
length scales near jc. As a result, the appropriate vari- 
able for a finite system near jc is tV-/'^ . Here we would 
like to emphasize that the function Ol appearing above 
depends on the lattice geometry as well as the boundary 
condition. For example, for the Ising model on a 2-D 
periodic square lattice, there exists a so-called shift ex- 
ponent A which is related to the deviation of the critical 
point on a finite volume from the bulk value. For the 
observables considered in this study, the most relevant 
FSS ansatz is given by ©^(t) = {I + cL-'^)goitL^^''), 
where the confiuent correction exponent co, which arises 
due to the inhomogeneous part of the free energy as well 
as the nonlinearity of the scaling field, is included ex- 
plicitly. One should be aware that above expression of 
FSS ansatz is an asymptotical one which is valid only for 
large L and close to jc- However to present the main 
results of this study, we find it is sufficient to employ 
the FSS ansatz introduced above for the data analysis. 
As mentioned earlier, when one approaches the critical 
point ( J'/ J)c for a second order phase transition, the ob- 
servables psi3A^^/^/2 and Q2 should satisfy the following 
finite-size scaling formula 

OL{t) = (l + cL-^)go{tL^^n 

= (1 + ci-") [50 + tL'^^gi + {tL'^n^92 

+ (tL'/n'93 + ■■■)], (5) 

where O stands for either psi3N^^^ /2 or Q2, t = {j — 
jc)/jc with j = J'/J and v and uj arc the correspond- 
ing critical exponents introduced before. Further, in 
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eq. ©, L = or L = N'^/'^ depending on whether 

the observable O is referring to PsiSN^^^ /2 or Q2. Fi- 
nally (7 is a smooth function of the variable tL'^ and 
50 + tL^^-'gi + {tL^'''?92 + {tL^'^faz + • • • is the Tay- 
lor expansion of the analytic function g. From eq. (O, 
one concludes that if the transition is second order, then 
the curves of different L should intersect at the critical 
point for large L. In the following, we will employ the 
finite-size scaling formula eq. ([5]) for various observables 
to calculate the critical point as well as the critical ex- 
ponent V. After determining the regime where the tran- 
sition occurs, we have performed substantial simulations 
with J' /J ranging from J'/J = 1.715 to J'/J = 1.755 
for N = 10^, 12^, 142, ... , 96^, where N is the total 
number of spins. We use sufficiently large (3 so that all 
the observables measured in our simulations take their 
zcro-tcmperature values. Figure [3] shows the data points 
for psiL and Q2 obtained from the simulations. From the 
figure, one observes that for both psiL and Q2 the curves 
of different L indeed have a tendency to intersect at one 
point for large L, which in turn is an indications for a 
second order phase transition. As a first step toward sys- 
tematically analyzing the data, we would like to have a 
better estimate of the critical point (J'/J)c. This can be 
obtained by investigating the crossing points from either 
PsiL or Q2 at system sizes L and 2L. Figured shows such 
crossing points as functions oil/ L. Indeed one sees both 
crossing points from psiL and Q2 approach a common 
J'/J with increasing L. By using the analysis suggested 
in (iGj . we find that {J'/J)c obtained from the upper 
(lower) crossing points in figure 3] is within the range 
iJ'/J)c e [1.7345,1.7365] {{J'/J)c & [1.7355,1.7363]) 
which agrees with the roughly estimated {J' /J)c from 
Xs- Next, we turn to determine the critical exponent v 
by employing the finite-size scaling ansatz for the observ- 
ables psiL and Q2- First of all, let us focus on Q2- A fit 
of Q2 to eq. (P with N > 12^ leads to (J'/J) = 1.7358(2) 
and ly = 0.691(2). The result is shown in figure[5l In the 
fits, we carefully make sure that v and (J'/J)c obtained 
from the fits are stable, namely v and [J'/J)c are con- 
sistent with respect to the elimination of smaller N. For 
example, from the fit using data points of = 24^ to 
N = 96^, we find (J7J)c = 1.7357(2) and ly = 0.694(2). 
Both of them agree with the results obtained from the fit 
using most of the available data points. With N > 26^, 
the fits are not stable. This is due to the fact that for 
large N, the data points of Q2 are not accurate enough to 
perform the fits including the term cL~^ in eq. jS]), where 
the corrections become very small. Therefore we use an- 
other approach, namely we use the data with sufficiently 
large A^ so that one can safely ignore the correction cL~'^ . 
We also make sure that the results are stable like what 
we did before. The validity of this strategy is justified 
as follows. We generate data points from the expression 
(1 -f- 0.5i-'')(l + {tL^I") + QAltL^'^f + 0.01(tLi/'=)3) for 
various values of L and t with 6 = 1.5 and c = 0.8. 
By analyzing this set of data, we do observe the con- 
vergence of c to c = 0.8 as one eliminates more data 



points with smaller value of L. Surprisingly, with this 
strategy we find that for a wide range of A^ we obtain 
consistent results with our previous analysis including 
the correction term cL^". For instance, using the data 
with A^ > 44^, we find v = 0.693(4) which agrees with 

V = 0.691(2) obtained earlier (figure [6]). At this stage, 
one naturally would conclude that we find the same un- 
conventional critical behavior as that in [l[. However, 
when we use larger A^, we begin to observe that v con- 
verges to the expected 0(3) value v = 0.7112(5). For 
example, fitting the data from A^ = 60^ to A^ = 96^ to 
eq. (O without the term cL~'^ leads to = 0.705(5), 
which is compatible with the expected 0(3) value (fig- 
ure [7|). By applying the same analysis to the observable 
PsiL, we arrive at similar results. For example, although 
with a slightly large x^/d.o.f. ~ 5.1, by fixing v = 1.7356 
which is the expected [J' / J)c calculated from Q2, wc find 

V = 0.694(4) for A^ > 48^ from the fit including the term 
cL^'^ in eq. ([5]). We attribute the poor quality of the 
fit to the correction to the critical point which are not 
considered here. Indeed with a fixed (J'/J)c = 1.7351 
which is not within the the statistical error of (J'/J)c 
obtained from Q2, we arrive a,t v = 0.707(3) and a better 
X^/d.o.f. = 3.7. Notice this v is even compatible with 
the known 0(3) value. With a fixed (J'/jjc = 1.7356, 
only from A^ > 64^ we begin to obtain v consistent with 
the expected 0(3) value v = 0.7112(5) : by fitting data 
points with A^ > 68^ to eq. (O with the term cL^'^ , we 
obtain ly = 0.707(6) (figured]). 

Table 1 summarizes the results of all the fits mentioned 
earlier. We would like to point out that the statistical er- 
rors shown in this study are obtained by binning. Hence 
the error bars presented here might be less accurately 
determined compared to those calculated with more so- 
phisticated methods. Indeed from the x^/d.o.f. shown 
in table 1, one sees that while the errors for Q2 seem 
to be overestimated, the errors for psiL are likely un- 
derestimated. Nevertheless, from the x^/d.o.f. presented 
in table 1, the results we have concluded should remain 
valid even an optimization procedure is applied for the 
data analysis. Notice the exponent oj for the last 2 rows 
in table 1 are not shown explicitly since the correspond- 
ing error bars are of the same magnitude as uj themselves. 
This is likely due to the fact that the data is not precise 
enough to determine tu unambiguously. Further, the uj 
we obtain from the fits are not consistent with the theo- 
retical value. For example, the uj in the last row of table 
1 is 0.38 which is significantly smaller than the theoreti- 
cal value OJ ~ 0.78. However since the value of cu will be 
affected by higher order operators to some extent as well 
as there might be a strong correlation between c and uj 
in formula ([5]), the uj presented here should be treated as 
a effective one. Actually we find that both x^/d.o.f. and 
ly from the fit of the last row in table 1 are very stable 
with a fixed value of uj including uj = 0.78. Therefore 
the result oi ly = 0.707(6) shown in table 1 should be 
trustable. 
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FIG. 4: Crossing points of paiL and Q2 obtained from the 
corresponding observables at system sizes L and 2L. The 
lines are added to guide the eye. 



TABLE I: Results of the fits mentioned in the text. A'^min and 
Aniax refer to the smallest and largest number of spins in the 
fit, respectively. 



IV. DISCUSSION AND CONCLUSION 

In this note, we have performed large scale Monte 
Carlo simulations to study the critical behavior of the 
spin- 1/2 Hcisenberg model with a spatially staggered 
anisotropy on the honeycomb lattice. In particular, we 
have demonstrated the subtleties of determining the crit- 
ical exponent v from investigating the finite-size scaling 
behavior of the relevant observables. From what we have 
found, we conclude that for the second order phase transi- 
tion considered here, the large volume data is essential to 
calculate the corresponding critical exponent v correctly. 
With a detailed data analysis, we find that the critical 
exponent v is given by = 0.707(6) which is cosistcnt 
with the most accurate 0(3) value v = 0.7112(5) ob- 
tained in [l3l ■ By including the sub- leading correction in 
the finite-size scaling ansatz and using most of the avail- 
able data points, we can obtain a value of v consistent 
with that calculated in [ij . Concerning the finite volume 
effect, we believe the v obtained from using only large 
N is more reliable. Our results do not necessarily imply 
that the critical exponent v obtained in [l[ should not be 
trusted. However from what we have found in this study, 
it would be desirable to go beyond the volumes used in 
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FIG. 5: Fit of most of the numerical data points for Q2 to 
their finite-size scaling ansatz including the correction term 
cL^'^ . While the open circles are the Monte Carlo data, the 
solid lines are obtained by using the results of the fit. Some 
data points are omitted for better visibility. 




FIG. 6: Fit of Q2 to its finite-size scaling ansatz given by 
eq. © without taking the correction cL~" into account. The 
number of spins N used in the fit ranges from 44^^ to 96^^. 
While the open circles are the Monte Carlo data, the solid 
lines are obtained by using the results of the fit. 



[l[ to see whether the unconventional behavior observed 
in [l| will continue to hold with the addition of larger 
volumes or not. 
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V = 0.705(5) 
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FIG. 7: Fit of the numerical data for Q2 to tlieir finite-size 
scaling ansatz without the correction term cL~^ in eq. ([Sjl. 
The number of spins N used in the fit ranges from A'' = 60^ 
toN ^ 96^ While the open circles are the Monte Carlo data, 
the solid lines are obtained by using the results of the fit. 
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FIG. 8: Fit of the numerical data for psiL to their finite-size 
scaling ansatz with the correction term cL~^ . While the open 
circles are the Monte Carlo data, the solid lines are obtained 
by using the results of the fit. 
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